El presente tema presenta algunos métodos de análisis en R.
Los temas cubiertos serán:
Al finalizar la sesión usted tendrá una noción de una gama de análisis que se pueden realizar en R, esperando que incentive a la persona a seguir buscando otras formas, métodos y técnicas de análisis de la información.
El muestreo es el proceso de seleccionar un conjunto de individuos de una población con el fin de estudiarlos y poder caracterizar el total de la población. Dada la imposibilidad de estudiar o análisis a una población entera (censar), debemos pasar el una muestreo, o por el método de muestreo.
El muestreo se constituye de dos partes:
No sé cubrirá el tamaño de muestre (n), y sus distintos diseños (aleatorio, sistemático, estratificado, conglomerado o complejo, etc) pero veremos dos formas de seleccionar casos en R.
Podemos tomar una muestra de tamaño n mediante la función “sample_n”. Tomemos una muestra de tamaño 50.
dim(iris)
## [1] 150 5
n <- iris %>% sample_n(50,replace=FALSE)
dim(n)
## [1] 50 5
También podemos seleccionar de forma aleatoria las observaciones pero de forma relativo o por una proporción.
Utilizamos la función “sample_frac”. Tomemos un 10% de la población.
dim(iris)
## [1] 150 5
n2 <- iris %>% sample_frac(0.1,replace=FALSE)
dim(n2)
## [1] 15 5
R posee dos librerías las siguientes librerías para tratar los temas de muestreo:
Mientras que sampling trabaja diseños aleatorios irrestrictos, sistemáticos y jerárquicos, survey se utiliza para diseños complejos o por conglomerados.
Abordamos dos tipos de enfoque:
Antes:
Una definición de la web explica:
Aprendizaje no supervisado es un método donde un modelo se ajusta a las observaciones. Se distingue del Aprendizaje supervisado por el hecho de que no hay un conocimiento a priori. En el aprendizaje no supervisado, un conjunto de datos de objetos de entrada es tratado.
En mis palabras, diría:
Empecemos primero los análisis de aprendizaje no supervisado.
El análisis por componente principales hace parte de un grupo de análisis descriptivos multidimensionales llamados métodos factoriales.
Son métodos descriptivos, no se apoyan en modelos probabilísticos sino más bien de un modelo geométrico para mejorar la representación multidimencional.
A partir de una matriz rectangular de datos con p variables cuantitativas y n unidades, el análisis por componentes principales propone diversas representación geométricas para el entendimiento de los individuos y las variables.
Lo que se busca es ver si existe una estructura, no conocida a priori, para el conjunto de casos y variables, y así mejorar la interpretación.
Como todo método descriptivo, llevar a cabo un PCA no es un fin en sí. El PCA sirve para conocer mejor los datos, detectar valores sospechosos, y ayuda a formular hipótesis que se deben estudiar luego mediante modelos predictivos.
Para la aplicación del PCA en R, se pueden utilizar diversas librerías y funciones. Se exponen algunos ejemplos:
Veamos la estrucuta de las salidas del análisis por PCA.
Veamos el porcentaje de variancia explicada por componente
Proporció de la varianza explicada
Proporció de la varianza explicada acumulada
El PCA busca mediante proyección explicar a:
Proyección de los individuos
Proyección de las variables
Proyección de individuos y variables
En el AF la atención se centra principalmente en las variables estadísticas y no en los individuos; es más un método de análisis multivariante que un método de análisis multidimensional.
El Análisis Factorial es un nombre genérico que se da a una clase de métodos estadísticos multivariantes cuyo propósito principal es definir la estructura subyacente en una matriz de datos. Generalmente hablando, aborda el problema de cómo analizar la estructura de las interrelaciones (correlaciones) entre un gran número de variables.
Al igual que el PCA, en el AF se busca reducir un número p a k variables. Sin embargo, el AF posee ciertas características que lo hacen el análisis ideal para conocer las estructura subyacente que se desea estudiar. Se dice que el AF es más elaborado, y se buscan conclusiones más contundentes, el PCA busca más aproximaciones a nivel gráfico del conjunto de datos. A esa reducción les llamamos Factores.
De forma matemática, el PCA es un AF pero más simple: sin extracción de vector propio, sin rotaciones, y sin transformaciones, además del gran supuesto practicamente irreal de la ortogonalidad de las variables o componentes principales en el PCA.
Desde un punto de vista práctico, el PCA es observacional, y el AF es un enfoque para un modelo explicito (estructura subyacente).
Mientras que el PCA busca el análisis de la variancia, el AF además busca el análisis de la covariancia y correlación entre las variables.
El objetivo típico en el análisis factorial es identificar las variables que están relacionadas entre sí, y separarlas de otras (una forma de agrupación de las variable). Los tipos de rotaciones ayudan a visualizar mejor este proceso.
Existen dos grandes tipos de modalidaes para el FA:
Veamos unos usos del FA exploratorio.
La forma más popular de conocer la cantidad óptima de factores es mediante el gráfico de sedimentación. Lo usual es ver dónde es que se lleva a cabo la caída más importante, o se dibuja un “codo”. A veces no es nada fácil el saber lo anterior. Personalmente considero que ante la poca claridad, hay que llevar a cabo un análisis interactivo del FA.
¿Cuántos factores?
R posee ciertos criterios que permiten “guiar” al analista en la elección del número de factor a elegir.
## Parallel analysis suggests that the number of factors = 3 and the number of components = NA
## noc naf nparallel nkaiser
## 1 5 1 5 5
Entonces, ¿cuántos componentes debemos elegir?
fa.1 <- fa(r=wine.1, nfactors = 3, rotate = "varimax")
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
fa.0 <- fa(r=personality, nfactors = 5, rotate = "varimax")
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
res1a <- factanal(personality, factors = 5, rotation = "varimax", na.action = na.omit)
# Calculadon la singularidad y la comunalidad
loadings_distant = res1a$loadings[1,]
communality_distant = sum(loadings_distant^2); communality_distant
## [1] 0.4554525
uniqueness_distant = 1-communality_distant; uniqueness_distant
## [1] 0.5445475
Veamos el peso de cada variable a los Factores
Factor 1
Factor 2
Factor 3
Según la cantidad de factores, la función fa.diagram nos puede indicar cuál es la estructura factorial del análisis.
LLevado a cabo el proceso de elegir la rotación y luego el número de componentes, podemos verificar que desde un inicio hasta la solución del problema, la estructura factorial del FA posee una forma de explicar los facores en un plano 2D. Veamos esto
Inicio: todas las variables
Intermedio: aplicación del FA y el acomodo de las variables
Final: estructura final del FA
Creamos la estructura factorial mediante el análisis de las cargas factoriales
shy = rowMeans(cbind(personality$distant, personality$shy, personality$withdrw, personality$quiet))
outgoing = rowMeans(cbind(personality$talkatv, personality$outgoin, personality$sociabl))
hardworking = rowMeans(cbind(personality$hardwrk, personality$persevr, personality$discipl))
friendly = rowMeans(cbind(personality$friendl, personality$kind, personality$coopera, personality$agreebl, personality$approvn, personality$sociabl))
anxious = rowMeans(cbind(personality$tense, personality$anxious, personality$worryin))
# Combinando los factores y creando una nueva estructura de datos
combined_data = cbind(shy,outgoing,hardworking,friendly,anxious)
combined_data = as.data.frame(combined_data)
Puede que acá se estimen diversos modelos de FA, con diversas rotaciones, con tal de ver una mejorar separación, y hasta aglomeración de los factores.
Es una técnica descriptiva o exploratoria cuyo objetivo es resumir una gran cantidad de datos en un número reducido de dimensiones, con la menor pérdida de información posible. En esta línea, su objetivo es similar al de los métodos factoriales, salvo que en el caso del análisis de correspondencias el método se aplica sobre variables categóricas (nominales y ordinales).
El análisis de correspondencias simples se utiliza a menudo en la representación de datos que se pueden presentar en forma de tablas de contingencia de dos variables nominales u ordinales. Otras utilizaciones implican el tratamiento de tablas de proximidad o distancia entre elementos, y tablas de preferencias.
El análisis de correspondencias consiste en resumir la información presente en las filas y columnas de manera que pueda proyectarse sobre un subespacio reducido, y representarse simultáneamente los puntos fila y los puntos columna, pudiéndose obtener conclusiones sobre relaciones entre las dos variables nominales u ordinales de origen.
La extensión del análisis de correspondencias simples al caso de varias variables nominales (tablas de contingencia multidimensionales) se denomina Análisis de Correspondencias Múltiples, y utiliza los mismos principios generales que la técnica anterior. En general se orienta a casos en los cuales una variable representa ítems o individuos y el resto son variables cualitativas u ordinales que representan cualidades.
housetasks
## Wife Alternating Husband Jointly
## Laundry 156 14 2 4
## Main_meal 124 20 5 4
## Dinner 77 11 7 13
## Breakfeast 82 36 15 7
## Tidying 53 11 1 57
## Dishes 32 24 4 53
## Shopping 33 23 9 55
## Official 12 46 23 15
## Driving 10 51 75 3
## Finances 13 13 21 66
## Insurance 8 1 53 77
## Repairs 0 3 160 2
## Holidays 0 1 6 153
res.ca <- CA(housetasks, graph = FALSE)
get_ca_row(res.ca)
## Correspondence Analysis - Results for rows
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the rows"
## 2 "$cos2" "Cos2 for the rows"
## 3 "$contrib" "contributions of the rows"
## 4 "$inertia" "Inertia of the rows"
get_ca_col(res.ca)
## Correspondence Analysis - Results for columns
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the columns"
## 2 "$cos2" "Cos2 for the columns"
## 3 "$contrib" "contributions of the columns"
## 4 "$inertia" "Inertia of the columns"
#--- para las filas:
fviz_contrib(res.ca, choice = "row", axes = 1)
#--- para las columnas:
fviz_contrib(res.ca, choice = "col", axes = 1)
#--- para las variables en las filas:
fviz_ca_row(res.ca, repel = TRUE)
#--- para las variables en las columnas:
fviz_ca_col(res.ca)
#--- y la visualización conjunta:
fviz_ca_biplot(res.ca, repel = TRUE)
El Análisis por Conglomerados (Cluster Analysis), es una técnica estadística multivariante que agrupa elementos (sujetos) tratando de lograr la máxima homogeneidad en cada grupo y la mayor diferencia entre los grupos.
Es un método basado en criterios geométricos y se utiliza fundamentalmente como una técnica exploratoria, descriptiva pero no explicativa. Se fundamenta en teorías de optimizaciones numéricas, y no en optimizaciones algebraicas.
Las soluciones no son únicas, en la medida en que la pertenencia al conglomerado para cualquier número de soluciones depende de muchos elementos del procedimiento elegido.
La solución del Cluster depende totalmente de las variables utilizadas, la adición o destrucción de variables relevantes puede tener un impacto substancial sobre la solución resultante.
Los algoritmos de formación de conglomerados se agrupan en diversas categorías, y existen muchas modalidades. El presente capítulo estudia los modalidades de clusters jerárquicos y los de k-medias.
El objetivo del análisis de cluster es, a partir de “𝑛” individuos (𝑛 grupos iniciales) buscar agrupar el conjunto de individios en dos o más grupos basados en la similitud de esos objetos según una serie especificada de características (variables).
Finalmente, siempre hay que tomar en cuenta los mismos tres aspectos que podrían influir fuertemente en el análisis de los datos:
Veremos algunos ejemplos de los algoritmos presentados.
En el entorno de programación R existen múltiples paquetes que implementan algoritmos de clustering y funciones para visualizar sus resultados. En este documento se emplean los siguientes:
stats: contiene las funciones:
Para evitar problemas de la magnitud de las variables, recordar SIEMPRE escalar las variables…
Veamos la estructura de algunas de funciones anteriores:
La estructura es :
#scale(x, center = TRUE, scale = TRUE)
La estructura es :
#kMeans(x, centers, iter.max=10, num.seeds=10)
La estructura es :
# fviz_cluster(object, data = NULL, choose.vars = NULL, stand = TRUE,
# axes = c(1, 2), geom = c("point", "text"), repel = FALSE,
# show.clust.cent = TRUE, ellipse = TRUE, ellipse.type = "convex",
# ellipse.level = 0.95, ellipse.alpha = 0.2, shape = NULL,
# pointsize = 1.5, labelsize = 12, main = "Cluster plot", xlab = NULL,
# ylab = NULL, outlier.color = "black", outlier.shape = 19,
# ggtheme = theme_grey(), ...)
La estructura es :
# hclust(d, method = "complete", members = NULL)
# S3 method for hclust
# plot(x, labels = NULL, hang = 0.1, check = TRUE,
# axes = TRUE, frame.plot = FALSE, ann = TRUE,
# main = "Cluster Dendrogram",
# sub = NULL, xlab = NULL, ylab = "Height", …)
El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Se pretende estudiar si existe una agrupación subyacente de los estados empleando K-means-clustering.
El paquete factoextra creado contiene funciones que facilitan en gran medida la visualización y evaluación de los resultados de clustering. Si se emplea K-means-clustering con distancia euclídea hay que asegurarse de que las variables empleadas son de tipo continuo, ya que trabaja con la media de cada una de ellas.
Veamos una representación de los datos, para 4 clusters (K=4)
datos <- scale(USArrests)
set.seed(123)
km_clusters <- kmeans(x = datos, centers = 4, nstart = 50)
# Las funciones del paquete factoextra emplean el nombre de las filas del
# dataframe que contiene los datos como identificador de las observaciones.
# Esto permite añadir labels a los gráficos.
fviz_cluster(object = km_clusters, data = datos, show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE) +
labs(title = "Resultados clustering K-means") +
theme_bw() +
theme(legend.position = "none")
Hierarchical clustering es una alternativa a los métodos de partitioning clustering que no requiere que se pre-especifique el número de clusters. Los métodos que engloba el hierarchical clustering se subdividen en dos tipos dependiendo de la estrategia seguida para crear los grupos:
Agglomerative clustering (bottom-up): el agrupamiento se inicia en la base del árbol, donde cada observación forma un cluster individual. Los clusters se van combinado a medida que la estructura crece hasta converger en una única “rama” central.
Divisive clustering (top-down): es la estrategia opuesta al agglomerative clustering, se inicia con todas las observaciones contenidas en un mismo cluster y se suceden divisiones hasta que cada observación forma un cluster individual.
En ambos casos, los resultados pueden representarse de forma muy intuitiva en una estructura de árbol llamada dendrograma.
En este se debe antes transformar los datos a distancias
# Crear las distancia y definir el método
dd <- dist(scale(USArrests), method = "euclidean")
hc <- hclust(dd, method = "ward.D2")
Veamos diversas modalidaes
Mediante el plot.hclust(), el más sencillo
plot(hc)
# Put the labels at the same height: hang = -1
plot(hc, hang = -1, cex = 0.6)
Mediante el plot.dendogram(),
# Convert hclust into a dendrogram and plot
hcd <- as.dendrogram(hc)
# Default plot
plot(hcd, type = "rectangle", ylab = "Height")
# Zoom in to the first dendrogram
plot(hcd, xlim = c(1, 20), ylim = c(1,8))
# Triangle plot
plot(hcd, type = "triangle", ylab = "Height")
Mediante el Phylogenetic trees
Este presenta diversas formas de representar los clusters
# Default plot
plot(as.phylo(hc), cex = 0.6, label.offset = 0.5)
# Cladogram
plot(as.phylo(hc), type = "cladogram", cex = 0.6,
label.offset = 0.5)
# Unrooted
plot(as.phylo(hc), type = "unrooted", cex = 0.6,
no.margin = TRUE)
# Fan
plot(as.phylo(hc), type = "fan")
# Radial
plot(as.phylo(hc), type = "radial")
# Change the appearance
# change edge and label (tip)
plot(as.phylo(hc), type = "cladogram", cex = 0.6,
edge.color = "steelblue", edge.width = 2, edge.lty = 2,
tip.color = "steelblue")
Ahora pasemos a los análisis supervisados.
En el aprendizaje supervisado una variable trata de ser explicada por las demás (relación asimétrica, variables dependientes e independientes).
La clasificación se centra en identificar a cuál de un conjunto de categorías (subpoblaciones) pertenece una nueva observación, sobre la base de un conjunto de datos de formación que contiene observaciones (o instancias) cuya categoría de miembros es conocida.
En el análisis de aprendizaje automático, para datos transversales, existen dos grandes categorias: la clasificación y la predicción.
Por ahora veremos la clasificación, luego en la regresión veremos la predicción, y volveremos a ver la clasifiación cuando vemos la regresión dicotómica.
El análisis discriminante es una técnica multivariante cuya finalidad es describir las diferencias significativas entre 𝑘 grupos de objetos (𝑘 > 1) sobre los que se observan (variable de clasificación).
En caso de que estas diferencias existan, intentará explicar en qué sentido se dan y proporcionar procedimientos de asignación sistemática de nuevas observaciones con grupo desconocido a uno de los grupos analizados, utilizando para ello sus valores en las𝑘variables clasificadoras.
Podemos ver este procedimiento como un modelo de predicción de una variable respuesta categórica (variable grupo) a partir de 𝑘 variables explicativas generalmente continuas (variables clasificatorias).
El Análisis Discriminante consiste en buscar nuevas variables (variables discriminantes), que separan lo mejor posible los grupos de proyección en las 𝑘 variables, de las 𝑛 observaciones, para poder predecir nuevas observaciones i.
En la clasificación se llevan a cado dos procesos fundamentales: construir la regla o ecuación de clasificación, y verificar la calidad del proceso, para luego introducir nuevos casos..
En la construcción de la regla de decisión, veremos un forma lineal y otro cuadrática.
Utilizaremos dos funciones lda() y qda() de la librería MASS.
La estructura es :
lda(x, …)
S3 method for formula
lda(formula, data, …, subset, na.action)
S3 method for default
lda(x, grouping, prior = proportions, tol = 1.0e-4,
method, CV = FALSE, nu, …)
S3 method for data.frame
lda(x, …)
S3 method for matrix
lda(x, grouping, …, subset, na.action)
Ver el enlace: https://www.rdocumentation.org/packages/MASS/versions/7.3-51.1/topics/lda
La estructura es :
qda(x, …)
S3 method for formula
qda(formula, data, …, subset, na.action)
S3 method for default
qda(x, grouping, prior = proportions,
method, CV = FALSE, nu, …)
S3 method for data.frame
qda(x, …)
S3 method for matrix
qda(x, grouping, …, subset, na.action)
Ver el enlace: https://www.rdocumentation.org/packages/MASS/versions/7.3-51.1/topics/qda
Pera evaluar la regla de decisión, se suele particionar el archivo de datos en entramiento y validación. Esto es un procedimiento típico de los análisis supervisados, y también en la minería de datos.
training_sample <- sample(c(TRUE, FALSE), nrow(iris), replace = T, prob = c(0.6,0.4))
train <- iris[training_sample, ]
test <- iris[!training_sample, ]
Dimensiones del set de entrenamiento
## [1] 89 5
Dimensiones del set de validación
## [1] 61 5
lda.iris <- lda(Species ~ ., train)
lda.iris
## Call:
## lda(Species ~ ., data = train)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3258427 0.3483146 0.3258427
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.051724 3.403448 1.462069 0.2344828
## versicolor 5.974194 2.783871 4.290323 1.3709677
## virginica 6.541379 2.979310 5.517241 1.9965517
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.8752702 -0.72840217
## Sepal.Width 1.4177443 3.31503628
## Petal.Length -2.0953556 0.07481247
## Petal.Width -3.0258487 1.42440523
##
## Proportion of trace:
## LD1 LD2
## 0.9911 0.0089
names(lda.iris)
## [1] "prior" "counts" "means" "scaling" "lev" "svd" "N"
## [8] "call" "terms" "xlevels"
Enfoquemos en la salida, la parte de los coeficientes (podemos verlo vediante el “scaling”)
¿Cuál sería entonces las ecuaciones de ambas reglas?
LD1 ? LD2 ?
qda.iris <- qda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, train)
qda.iris
## Call:
## qda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
## data = train)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3258427 0.3483146 0.3258427
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.051724 3.403448 1.462069 0.2344828
## versicolor 5.974194 2.783871 4.290323 1.3709677
## virginica 6.541379 2.979310 5.517241 1.9965517
Analisando la función qda, no se programó la ecuación del qda…
Una forma de ver la efectividad de las reglas de separarción en poner en un plano 2D los casos con sus respectivas reglas si hubiese 2.
plot(lda.iris, col = as.integer(train$Species))
También se puede corroborar mediante el análisis de la posición y variabilidad por categoría de clasificación, en cada regla. El caso de la RD1:
lda.iris <- lda(Species ~ ., train)
lda.iris
## Call:
## lda(Species ~ ., data = train)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3258427 0.3483146 0.3258427
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.051724 3.403448 1.462069 0.2344828
## versicolor 5.974194 2.783871 4.290323 1.3709677
## virginica 6.541379 2.979310 5.517241 1.9965517
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.8752702 -0.72840217
## Sepal.Width 1.4177443 3.31503628
## Petal.Length -2.0953556 0.07481247
## Petal.Width -3.0258487 1.42440523
##
## Proportion of trace:
## LD1 LD2
## 0.9911 0.0089
plot(lda.iris, dimen = 1, type = "b")
Tengo un error para la 2nd dimensión… creo que debe ser porque explica muy poco.
#plot(lda.iris, dimen = 2, type = "both")
Otra forma de ver la efectividad de las reglas es mediante la proyección de la regla y sus inviduos para ciertas variables.
Veamos que tanto separa o discrimina la primera regla
Y para la 2nda regla:
¿Qué podemos concluir?
El uso de la función partimat del paquete klaR proporciona una forma alternativa de trazar las funciones discriminantes lineales. partimat emite una serie de gráficos para cada combinación de dos variables. Piense en cada gráfico como una vista diferente de los mismos datos. Las regiones coloreadas delinean cada área de clasificación. Se predice que cualquier observación que caiga dentro de una región será de una clase específica. Cada gráfico también incluye la tasa de error aparente para esa vista de los datos.
Partición lineal
Partición cuadrática
La función partimat() del paquete klar permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.
Partición lineal
Partición cuadrática
Podemos verificar la calidad del análisis mediante las proyeccioes de las variables y las predicciones.
Variables: podemos ver los resultados del modelo de la separación de los casos en los variables.
¿Qué podemos decir?
Predicciones.
Podemos para los valores de las reglas verificar si los individuos se están o no separando o aglomerando en cierto parte.
LD1: ¿qué podemos decir?
LD2: ¿qué podemos decir?
##
## setosa versicolor virginica
## setosa 29 0 0
## versicolor 0 30 1
## virginica 0 1 28
–> Recomendable tenerlo mejor en %
¿Qué podemos concluir del modelo para el entrenamiento?
##
## setosa versicolor virginica
## setosa 29 0 0
## versicolor 0 30 1
## virginica 0 1 28
–> Recomendable tenerlo mejor en %
¿Qué podemos concluir del modelo para el entrenamiento? Y comparando LDA y QDA
##
## setosa versicolor virginica
## setosa 21 0 0
## versicolor 0 18 0
## virginica 0 1 21
–> Recomendable tenerlo mejor en %
¿Qué podemos concluir del modelo para el entrenamiento?
##
## setosa versicolor virginica
## setosa 21 0 0
## versicolor 0 18 0
## virginica 0 1 21
–> Recomendable tenerlo mejor en %
¿Qué podemos concluir del modelo para el entrenamiento? Y comparando LDA y QDA
El objetivo es mostrar los principales comandos en R para generar:
De igual forma se verán ciertas etapas de análisis que componen a estas variantes de la regresión.
El dataset Boston del paquete MASS recoge la mediana del valor de la vivienda en 506 áreas residenciales de Boston. Junto con el precio, se han registrado 13 variables adicionales.
tail(Boston)
## crim zn indus chas nox rm age dis rad tax ptratio black lstat
## 501 0.22438 0 9.69 0 0.585 6.027 79.7 2.4982 6 391 19.2 396.90 14.33
## 502 0.06263 0 11.93 0 0.573 6.593 69.1 2.4786 1 273 21.0 391.99 9.67
## 503 0.04527 0 11.93 0 0.573 6.120 76.7 2.2875 1 273 21.0 396.90 9.08
## 504 0.06076 0 11.93 0 0.573 6.976 91.0 2.1675 1 273 21.0 396.90 5.64
## 505 0.10959 0 11.93 0 0.573 6.794 89.3 2.3889 1 273 21.0 393.45 6.48
## 506 0.04741 0 11.93 0 0.573 6.030 80.8 2.5050 1 273 21.0 396.90 7.88
## medv
## 501 16.8
## 502 22.4
## 503 20.6
## 504 23.9
## 505 22.0
## 506 11.9
summary(Boston)
## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00
Veamos las distribuciones por variable:
multi.hist(x = Boston[,1:3], dcol = c("blue","red"), dlty = c("dotted", "solid"),
main = "")
multi.hist(x = Boston[,5:9], dcol = c("blue","red"), dlty = c("dotted", "solid"),
main = "")
multi.hist(x = Boston[,10:14], dcol = c("blue","red"),
dlty = c("dotted", "solid"), main = "")
vemos que los datos no poseen total normalidad en las variables.
Se pretende predecir el valor de la vivienda en función del porcentaje de pobreza de la población. Empleando la función lm() se genera un modelo de regresión lineal por mínimos cuadrados en el que la variable respuesta es medv y el predictor lstat.
modelo_simple <- lm(data = Boston,formula = medv ~ lstat)
modelo_simple
##
## Call:
## lm(formula = medv ~ lstat, data = Boston)
##
## Coefficients:
## (Intercept) lstat
## 34.55 -0.95
names(modelo_simple)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
summary(modelo_simple)
##
## Call:
## lm(formula = medv ~ lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.168 -3.990 -1.318 2.034 24.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.55384 0.56263 61.41 <2e-16 ***
## lstat -0.95005 0.03873 -24.53 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared: 0.5441, Adjusted R-squared: 0.5432
## F-statistic: 601.6 on 1 and 504 DF, p-value: < 2.2e-16
En la información devuelta por el summary se observa que el p-value del estadístico F es muy pequeño, indicando que al menos uno de los predictores del modelo está significativamente relacionado con la variable respuesta. ¿Cómo sería en palabras del presente ejemplo?
La creación de un modelo de regresión lineal simple suele acompañarse de una representación gráfica superponiendo las observaciones con el modelo. Además de ayudar a la interpretación, es el primer paso para identificar posibles violaciones de las condiciones de la regresión lineal.
attach(Boston)
plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
abline(modelo_simple, lwd = 3, col = "red")
La representación gráfica de las observaciones muestra que la relación entre ambas variables estudiadas no es del todo lineal, lo que apunta a que otro tipo de modelo podría explicar mejor la relación. Aun así la aproximación no es mala.
Una vez generado el modelo, es posible predecir el valor de la vivienda sabiendo el estatus de la población en la que se encuentra. Toda predicción tiene asociado un error y por lo tanto un intervalo. Es importante diferenciar entre dos tipos de intervalo:
Intervalo de confianza: Devuelve un intervalo para el valor promedio de todas las viviendas que se encuentren en una población con un determinado porcentaje de pobreza, supóngase lstat=10.
predict(object = modelo_simple, newdata = data.frame(lstat = c(10)),
interval = "confidence", level = 0.95)
## fit lwr upr
## 1 25.05335 24.47413 25.63256
Intervalo de predicción: Devuelve un intervalo para el valor esperado de una vivienda en particular que se encuentre en una población con un determinado porcentaje de pobreza.
predict(object = modelo_simple, newdata = data.frame(lstat = c(10)),
interval = "prediction", level = 0.95)
## fit lwr upr
## 1 25.05335 12.82763 37.27907
Una de las mejores formas de confirmar que las condiciones necesarias para un modelo de regresión lineal simple por mínimos cuadrados se cumplen es mediante el estudio de los residuos del modelo.
En R, los residuos se almacenan dentro del modelo bajo el nombre de residuals. R genera automáticamente los gráficos más típicos para la evaluación de los residuos de un modelo.
par(mfrow = c(1,2))
plot(modelo_simple)
par(mfrow = c(1,1))
Otra forma de identificar las observaciones que puedan ser outliers o puntos con alta influencia (leverage) es emplear las funciones rstudent() y hatvalues().
plot(x = modelo_simple$fitted.values, y = abs(rstudent(modelo_simple)),
main = "Absolute studentized residuals vs predicted values", pch = 20,
col = "grey30")
abline(h = 3, col = "red")
plot(hatvalues(modelo_simple), main = "Medición de leverage", pch = 20)
# Se añade una línea en el threshold de influencia acorde a la regla
# 2.5x((p+1)/n)
abline(h = 2.5*((dim(modelo_simple$model)[2]-1 + 1)/dim(modelo_simple$model)[1]),
col = "red")
En este caso muchos de los valores parecen posibles outliers o puntos con alta influencia porque los datos realmente no se distribuyen de forma lineal en los extremos. Deberíamos evaluar otros modelos y otras variables.
Se desea generar un modelo que permita explicar el precio de la vivienda de una población empleando para ello cualquiera de las variables disponibles en el dataset Boston y que resulten útiles en el modelo.
R permite crear un modelo con todas las variables incluidas en un data.frame de la siguiente forma:
modelo_multiple <- lm(formula = medv ~ ., data = Boston)
# También se pueden especificar una a una
summary(modelo_multiple)
##
## Call:
## lm(formula = medv ~ ., data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## crim -1.080e-01 3.286e-02 -3.287 0.001087 **
## zn 4.642e-02 1.373e-02 3.382 0.000778 ***
## indus 2.056e-02 6.150e-02 0.334 0.738288
## chas 2.687e+00 8.616e-01 3.118 0.001925 **
## nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
## rm 3.810e+00 4.179e-01 9.116 < 2e-16 ***
## age 6.922e-04 1.321e-02 0.052 0.958229
## dis -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
## rad 3.060e-01 6.635e-02 4.613 5.07e-06 ***
## tax -1.233e-02 3.760e-03 -3.280 0.001112 **
## ptratio -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
## black 9.312e-03 2.686e-03 3.467 0.000573 ***
## lstat -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
En el summary se puede observar que algunos predictores tienen p-values muy altos, sugiriendo que no contribuyen al modelo por lo que deben ser excluidos, por ejemplo age e indus. La exclusión de predictores basándose en p-values no es aconsejable, en su lugar se recomienda emplear métodos de best subset selection, stepwise selection (forward, backward e hybrid) o Shrinkage/regularization.
step(modelo_multiple, direction = "both", trace = 0)
##
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat, data = Boston)
##
## Coefficients:
## (Intercept) crim zn chas nox rm
## 36.341145 -0.108413 0.045845 2.718716 -17.376023 3.801579
## dis rad tax ptratio black lstat
## -1.492711 0.299608 -0.011778 -0.946525 0.009291 -0.522553
La selección de predictores empleando stepwise selection (hybrid/doble) ha identificado como mejor modelo el formado por los predictores crim, zn, chas, nox, rm, dis, rad, tax, ptratio, black, lstat.
modelo_multiple <- lm(formula = medv ~ crim + zn + chas + nox + rm + dis +
rad + tax + ptratio + black + lstat, data = Boston)
# También se pueden indicar todas las variables de un data.frame y exluir algunas
# modelo_multiple <- lm(formula = medv~. -age -indus, data = Boston)
summary(modelo_multiple)
##
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat, data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5984 -2.7386 -0.5046 1.7273 26.2373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.341145 5.067492 7.171 2.73e-12 ***
## crim -0.108413 0.032779 -3.307 0.001010 **
## zn 0.045845 0.013523 3.390 0.000754 ***
## chas 2.718716 0.854240 3.183 0.001551 **
## nox -17.376023 3.535243 -4.915 1.21e-06 ***
## rm 3.801579 0.406316 9.356 < 2e-16 ***
## dis -1.492711 0.185731 -8.037 6.84e-15 ***
## rad 0.299608 0.063402 4.726 3.00e-06 ***
## tax -0.011778 0.003372 -3.493 0.000521 ***
## ptratio -0.946525 0.129066 -7.334 9.24e-13 ***
## black 0.009291 0.002674 3.475 0.000557 ***
## lstat -0.522553 0.047424 -11.019 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7348
## F-statistic: 128.2 on 11 and 494 DF, p-value: < 2.2e-16
En los modelos de regresión lineal con múltiples predictores, además del estudio de los residuos vistos en el modelo simple, es necesario descartar colinealidad o multicolinealidad entre variables.
par(mfrow = c(1,2))
plot(modelo_multiple)
par(mfrow = c(1,1))
Para la colinealidad se recomienda calcular el coeficiente de correlación entre cada par de predictores incluidos en el modelo:
require(corrplot)
corrplot.mixed(corr = cor(Boston[,c("crim", "zn", "nox", "rm", "dis", "rad",
"tax", "ptratio", "black", "lstat", "medv")],
method = "pearson"))
El análisis muestra correlaciones muy altas entre los predictores rad y tax (positiva) y entre dis y nox (negativa).
attach(Boston)
## The following objects are masked from Boston (pos = 3):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, zn
par(mfrow = c(2,2))
plot(x = tax, y = rad, pch = 20)
plot(x = tax, y = nox, pch = 20)
plot(x = dis, y = nox, pch = 20)
plot(x = medv, y = rm, pch = 20)
par(mfrow = c(1,1))
Si la correlación es alta y por lo tanto las variables aportan información redundante, es recomendable analizar si el modelo mejora o no empeora excluyendo alguno de estos predictores.
Para el estudio de la multicolinealidad una de las medidas más utilizadas es el factor de inflación de varianza VIF. Puede calcularse mediante la función vif() del paquete car.
require(car)
## Loading required package: car
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:boot':
##
## logit
## The following object is masked from 'package:purrr':
##
## some
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
vif(modelo_multiple)
## crim zn chas nox rm dis rad tax
## 1.789704 2.239229 1.059819 3.778011 1.834806 3.443420 6.861126 7.272386
## ptratio black lstat
## 1.757681 1.341559 2.581984
Los indices VIF son bajos o moderados, valores entre 5 y 10 indican posibles problemas y valores mayores o iguales a 10 se consideran muy problemáticos.
De igual forma, podemos realizar la predicción para el modelo múltiple pero indicando el valor a predir para cada variable.
La Regresión Polinomial, aunque permite describir relaciones no lineales, se trata de un modelo lineal en el que se incorporan nuevos predictores elevando el valor de los ya existentes a diferentes potencias.
Cuando se intenta predecir el valor de la vivienda en función del estatus de la población, el modelo lineal generado no se ajusta del todo bien debido a que las observaciones muestran una relación entre ambas variables con cierta curvatura.
attach(Boston)
## The following objects are masked from Boston (pos = 5):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, zn
## The following objects are masked from Boston (pos = 6):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, zn
plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
abline(modelo_simple, lwd = 3, col = "red")
La curvatura descrita apunta a una posible relación cuadrática, por lo que un polinomio de segundo grado podría capturar mejor la relación entre las variables. En R se pueden generar modelos de regresión polinómica de diferentes formas:
Identificando cada elemento del polinomio: modelo_pol2 <- lm(formula = medv ~ lstat + I(lstat^2), data = Boston) El uso de I() es necesario ya que el símbolo ^ tiene otra función dentro de las formula de R.
Con la función poly(): lm(formula = medv ~ poly(lstat, 2), data = Boston)
modelo_pol2 <- lm(formula = medv ~ poly(lstat, 2), data = Boston)
summary(modelo_pol2)
##
## Call:
## lm(formula = medv ~ poly(lstat, 2), data = Boston)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2834 -3.8313 -0.5295 2.3095 25.4148
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.5328 0.2456 91.76 <2e-16 ***
## poly(lstat, 2)1 -152.4595 5.5237 -27.60 <2e-16 ***
## poly(lstat, 2)2 64.2272 5.5237 11.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.524 on 503 degrees of freedom
## Multiple R-squared: 0.6407, Adjusted R-squared: 0.6393
## F-statistic: 448.5 on 2 and 503 DF, p-value: < 2.2e-16
El p-value próximo a 0 del predictor cuadrático de lstat indica que contribuye a mejorar el modelo.
attach(Boston)
## The following objects are masked from Boston (pos = 3):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, zn
## The following objects are masked from Boston (pos = 6):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, zn
## The following objects are masked from Boston (pos = 7):
##
## age, black, chas, crim, dis, indus, lstat, medv, nox, ptratio, rad,
## rm, tax, zn
plot(x = lstat, y = medv, main = "medv vs lstat", pch = 20, col = "grey30")
points(lstat, fitted(modelo_pol2), col = 'red', pch = 20)
A la hora de comparar dos modelos se pueden evaluar sus R2. En este caso el modelo cuadrático es capaz de explicar un 64% de variabilidad frente al 54% del modelo lineal.
Dado que un polinomio de orden n siempre va a estar anidado a uno de orden n+1, se pueden comparar modelos polinómicos dentro un rango de grados haciendo comparaciones secuenciales.
anova(modelo_simple, modelo_pol2)
## Analysis of Variance Table
##
## Model 1: medv ~ lstat
## Model 2: medv ~ poly(lstat, 2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 504 19472
## 2 503 15347 1 4125.1 135.2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El p-value obtenido para el estadístico F confirma que el modelo cuadrático es superior.
Si en vez de una variable continua tenemos una variable categórica, y esa desea ser explicada por un conjunto de predictores o variables independienntes, podemos aplicar una regresión dicotómica.
Existen dos grandes tipos:
r: indicador de presencia de oxalato de calcio. 0: ausencia de oxalato en orina.1: presencia de oxalato en orina.gravity: gravedad específica de la orina.ph: pH de la orina.osmo: osmolaridad de la orina.cond: conductivicad de la orina.urea: concentración de urea en la orina.calc: concentración de calcio (milimoles por litro).datos <- read.csv("C:/Users/oscar/Desktop/R --- SAF/Tema 7/Orina.csv")
# Conversión de variable objetivo a factor
datos$r <- as.factor(datos$r)
# Imprimiendo datos
tail(datos)
## r gravity ph osmo cond urea calc
## 74 1 1.022 5.09 736 19.8 418 8.53
## 75 1 1.025 7.90 721 23.6 301 9.04
## 76 1 1.017 4.81 410 13.3 195 0.58
## 77 1 1.024 5.40 803 21.8 394 7.82
## 78 1 1.016 6.81 594 21.4 255 12.20
## 79 1 1.015 6.03 416 12.8 178 9.39
\[Y_i \sim\ B(p_i,\ n_i),\ para\ i=1,...,m\]
\[p(Y =1\ |\ X)=\frac{e^{(\beta_0+\beta_1x)}}{e^{(\beta_0+\beta_1x)}+1}\]
\[p(X) = \frac{e^{(\beta_0+\beta_1x)}}{e^{(\beta_0+\beta_1x)}+1}\\ p(e^{(\beta_0+\beta_1x)}+1) = e^{(\beta_0+\beta_1x)}\\ p \times e^{(\beta_0+\beta_1x)}\ +\ p = e^{(\beta_0+\beta_1x)}\\ p = e^{(\beta_0+\beta_1x)}\ -\ p \times e^{(\beta_0+\beta_1x)}\\ p = e^{(\beta_0+\beta_1x)}(1-p)\\ \frac{p}{1-p} = e^{(\beta_0+\beta_1x)}\]
\[ln(\frac{p}{1-p}) = \beta_0+\beta_1x\]
\[Y_i \sim\ B(p_i,\ n_i),\ para\ i=1,...,m\]
\[p(Y = 1|X)=\ \Phi(X^{T}\beta)\]
library(tidyr)
datos %>%
gather(key = "variable", value = "valor", -r) %>%
group_by(variable, r) %>%
summarise(Promedio = round(mean(valor, na.rm = TRUE), digits = 2),
`D. Estándar` = round(sd(valor, na.rm = TRUE), digits = 2),
Mínimo = round(min(valor, na.rm = TRUE), digits = 2),
Máximo = round(max(valor, na.rm = TRUE), digits = 2),
Q1 = round(quantile(valor, prob = 0.25, na.rm = TRUE), digits = 2),
Q2 = round(quantile(valor, prob = 0.5, na.rm = TRUE), digits = 2),
Q3 = round(quantile(valor, prob = 0.75, na.rm = TRUE), digits = 2)) %>%
rename(Oxalato = r, Variable = variable) %>%
datatable()
## `summarise()` regrouping output by 'variable' (override with `.groups` argument)
Análisis Exploratorio
library(broom)
tidy(apply(datos, MARGIN = 2, is.na)) %>%
gather(key = "variable", value = "valor") %>%
mutate(valor = as.numeric(valor)) %>%
group_by(variable) %>%
summarise(Total_NAs = sum(valor))
## Warning: 'tidy.logical' is deprecated.
## See help("Deprecated")
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## `summarise()` ungrouping output (override with `.groups` argument)
## Warning: `...` is not empty.
##
## We detected these problematic arguments:
## * `needs_dots`
##
## These dots only exist to allow future extensions and should be empty.
## Did you misspecify an argument?
## # A tibble: 1 x 2
## variable Total_NAs
## <chr> <dbl>
## 1 x 0
library(ggplot2)
colores <- c("dodgerblue", "gray40")
datos %>%
rename(Oxalato = r) %>%
gather(key = "variable", value = "valor", -Oxalato) %>%
ggplot(data = ., aes(x = valor, fill = Oxalato)) +
facet_wrap(~variable, scales = "free") +
geom_density(alpha = 0.7) +
scale_fill_manual(values = colores) +
labs(x = "", y = "") +
theme_light() +
theme(legend.position = "bottom")
## Warning: Removed 2 rows containing non-finite values (stat_density).
par(mfrow = c(2, 3))
cdplot(datos$r ~ datos$gravity, xlab = "Gravedad específica",
ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$ph, xlab = "pH",
ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$osmo, xlab = "Osmolaridad",
ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$cond, xlab = "Conductividad",
ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$urea, xlab = "Urea",
ylab = "Oxalato", col = colores)
cdplot(datos$r ~ datos$calc, xlab = "Calcio",
ylab = "Oxalato", col = colores)
datos %>%
rename(Oxalato = r) %>%
gather(key = "variable", value = "valor", -Oxalato) %>%
ggplot(data = ., aes(x = Oxalato, y = valor, fill = Oxalato)) +
facet_wrap(~variable, scales = "free") +
geom_boxplot() +
scale_fill_manual(values = colores) +
labs(x = "", y = "") +
theme_light() +
theme(legend.position = "bottom")
## Warning: Removed 2 rows containing non-finite values (stat_boxplot).
# Modelos Lineales Generalizados
mod_logit <- glm(r ~ ., data = datos, family = binomial(link = "logit"))
mod_probit <- glm(r ~ ., data = datos, family = binomial(link = "probit"))
# -------- Validación LOOCV (Manual)
out_logi_0.5 <- NULL
out_prob_0.5 <- NULL
for(i in 1:nrow(datos)){
out_logi_0.5[i] = predict(update(mod_logit, data = datos[-i, ]),
newdata = datos[i,], type = "response")
out_prob_0.5[i] = predict(update(mod_probit, data = datos[-i, ]),
newdata = datos[i,], type = "response")
}
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
# -------- Validación LOOCV (cv.glm())
## Función de coste con cutoff = 0.5
coste_0.5 <- function(r, pi = 0) mean(abs(r-pi)> 0.5)
## LOOCV
library(boot)
cv_error_logi_0.5 <- cv.glm(data = datos %>% filter(!is.na(osmo) & !is.na(cond)),
glmfit = mod_logit, cost = coste_0.5)
cv_error_prob_0.5 <- cv.glm(data = datos %>% filter(!is.na(osmo) & !is.na(cond)),
glmfit = mod_probit, cost = coste_0.5)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
cv.glm() requiere una función de coste para calcular el error. El objeto devuelto por la función del paquete boot almacena el error con el nombre delta; al restar 1 menos el error (delta) se obtendrá la precisión o accuracy del modelo, que es exactamente el mismo valor obtenido manualmente.
Función de coste: en el código hay una función de coste o pérdida para el límite igual a 0.5. Esta función puede ser expresada de la siguiente manera: \(error = \sum |r_i - p_i| > 0.5\). Donde \(r_i\) es el i-ésimo valor real y \(p_i\) es el i-ésimo valor predicho. En el siguiente código es posible evidenciar que se obtienen los mismos resultados de forma manual y con la función cv.glm().
Los modelos dicotómicos también sirven para la predicción, y según cierta condición, brindarán una probabilidad (entre 0-1) sobre un evento que así se consulte a través del modelo estimado.
El análisis de las series temporales juega un rol esencial tanto para conocer un determinado fenómeno como para pronosticar hacia el futuro.
Las series de tiempo son observaciones sobre un determinado fenómeno efectuadas en el tiempo, en lapsos ojala equivalente, o con intervalos de igual valor.
Ejemplos: exportaciones mensuales, ventas diarias de un producto, casos semanales de sida, temperatura promedio diaria, tasa anual de mortalidad, numero mensual de divorcios.
Una serie de tiempo posee dos características esenciales:
Los valores están ordenados o presentados en el tiempo.
Existe una dependencia o correlación entre los valores de la serie en el tiempo.
Un análisis por series temporales podría buscar:
Descripción de la serie Adecuar un modelo o técnica estocástica Pronostico en un periodo h de la serie.
El análisis de la serie debe preguntarse sobre el tipo de serie que se está analizando, los tipos de datos, y el período de referencia en la adecuación del mejor modelo estocástico.
De igual forma, dependiendo de la serie, el pronóstico debe considerar ciertas restricciones.
El análisis de una serie de tiempo de proceder, sin ser restrictivo, procede a groso modo de la siguiente forma:
Para faciliar las cosas, hablemos de la regresión en un contexto transversal y otro longitudinal.
En un contexto transversal, nuestro puntos o sujetos de estudio son las personas en un momento ti. En un contexto longitudinal, nuestro puntos o sujetos de estudio son los peridos de todo un rango de períodos t.
Mientras que la predicción busca saber el efecto de los coeficientes y realizar una intrapolación, en el pronóstico buscamos lo opuesto: extrapolamos y no le damos importancia al valor de los coeficientes.
Primero cargamos los datos y los transformados en ts (time series).
co2 <- read.csv("C:/Users/oscar/Desktop/R --- SAF/Tema 7/co2.csv", sep=";")
co2ts<-ts(co2$CO2, start = c(1959,1), frequency = 12)
La estructura de una serie ts:
print(co2ts)
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1959 333.13 332.09 331.10 329.14 327.36 327.29 328.23 329.55 330.62 331.40
## 1960 333.92 333.43 331.85 330.01 328.51 328.41 329.25 330.97 331.60 332.60
## 1961 334.68 334.17 332.96 330.80 328.98 328.57 330.20 331.58 332.67 333.17
## 1962 336.82 336.12 334.81 332.56 331.30 331.22 332.37 333.49 334.71 335.23
## 1963 337.95 338.00 336.37 334.47 332.46 332.29 333.76 334.80 336.00 336.63
## 1964 339.05 339.27 337.64 335.68 333.77 334.09 335.29 336.76 337.77 338.26
## 1965 341.47 341.31 339.41 337.74 336.07 336.07 337.22 338.38 339.32 340.41
## 1966 343.02 342.54 340.88 338.75 337.05 337.13 338.45 339.85 340.90 341.70
## 1967 344.28 343.42 342.02 339.97 337.84 338.00 339.20 340.63 341.41 342.68
## 1968 345.92 345.40 344.16 342.11 340.11 340.15 341.38 343.02 343.87 344.59
## 1969 347.38 346.78 344.96 342.71 340.86 341.13 342.84 344.32 344.88 345.62
## 1970 348.53 347.87 346.00 343.86 342.55 342.57 344.11 345.49 346.04 346.70
## 1971 349.93 349.26 347.44 345.55 344.21 343.67 345.09 346.27 347.33 347.82
## 1972 351.71 350.94 349.10 346.77 345.73
## Nov Dec
## 1959 331.87 333.18
## 1960 333.57 334.72
## 1961 334.86 336.07
## 1962 336.54 337.79
## 1963 337.93 338.95
## 1964 340.10 340.88
## 1965 341.69 342.51
## 1966 342.70 343.65
## 1967 343.04 345.27
## 1968 345.11 347.07
## 1969 347.23 347.62
## 1970 347.38 349.38
## 1971 349.29 350.91
## 1972
Veamos la serie como tal:
autoplot(co2ts, ts.colour = "blue", ts.linetype = "dashed")
¿Qué podemos decir? ¿Hay autocorrelación y estacionalidad en el tiempo?
Un gráfico de los autocorrelaciones totales nos puede ayudar:
autoplot(acf(co2ts, plot = FALSE))
## Warning: `mutate_()` is deprecated as of dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
Vemos un estudio de los componentes de una serie:
autoplot(stl(co2ts, s.window = "periodic"), ts.colour = "blue")
Vamos a estimar un total de 6 modelos ARIMA y luego comprar al mejor entre ellos:
arima1<- Arima(co2ts, order=c(0,1,2), seasonal=list(order=c(0,1,1),period=12))
arima2<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(2,1,0),period=12))
arima3<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(2,1,1),period=12))
arima4<- Arima(co2ts, order=c(1,1,1), seasonal=list(order=c(2,1,1),period=12))
arima5<- Arima(co2ts, order=c(1,1,2), seasonal=list(order=c(1,1,1),period=12))
arima6<- Arima(co2ts, order=c(0,1,1), seasonal=list(order=c(0,1,1),period=12))
arima7<- Arima(co2ts, order=c(1,1,0), seasonal=list(order=c(1,1,0),period=12))
Veamos las salidas del ARIMA(0,1,2)(0,1,1)
arima1
## Series: co2ts
## ARIMA(0,1,2)(0,1,1)[12]
##
## Coefficients:
## ma1 ma2 sma1
## -0.5126 0.0272 -0.7499
## s.e. 0.0813 0.0854 0.0920
##
## sigma^2 estimated as 0.1041: log likelihood=-38.86
## AIC=85.73 AICc=86.01 BIC=97.71
¿Qué nos dice la salida? ¿Cómo interpretamos los coeficientes?
Veamos cuál es el mejor modelo entre los 6 estimados. Utilizaremos el criterio de AIC y BIC:
AIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)
## df AIC
## arima1 4 85.72537
## arima2 4 93.70131
## arima3 7 90.55929
## arima4 6 89.02326
## arima5 6 88.62985
## arima6 3 83.82626
## arima7 3 106.19086
BIC(arima1,arima2,arima3,arima4,arima5,arima6,arima7)
## df BIC
## arima1 4 97.71422
## arima2 4 105.69016
## arima3 7 111.53978
## arima4 6 107.00654
## arima5 6 106.61312
## arima6 3 92.81789
## arima7 3 115.18250
Aquí se puede apreciar que los ajustes que mejor AIC y BIC presentan son aquellas que solo tienen componente de médias móviles y no tienen componente autorregresiva, siendo ARIMA(0,1,1)(0,1,1) el modelo que los test arrojan con un menor valor y por tanto con una mayor consideración.
Una vez estimados los modelos y elegido el mejor de ellos, en este caso ARIMA(0,1,1)(0,1,1), se procede a validarlo, esto es, ver la adecuidad del modelo seleccionado. Para ello en primer lugar se grafican los correlogramas de los residuos para comprobar que son ruido blanco:
ggtsdiag(arima6)
El modelo respecta los criterior de adecuación.
Finalizada la etapa de estimación, selección y validación de un modelo candidato, vamos a proyectar para los siguientes cuatro años siguientes, lo cuál son 48 meses, o un horizonte de proyección h=48.
forecast1<-forecast(arima6, level = c(95), h = 50)
autoplot(forecast1)
## Warning: `filter_()` is deprecated as of dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
En ete caso solo se detalló una estimación por un modelo ARIMA, pero para series univariadas se podría optar por: